# Bernd Beber, Michael J. Gilligan, Jenny Guardado, and Sabrina Karim
# R code to replicate analysis for "Peacekeeping, International Norms, and Transactional Sex in Monrovia, Liberia"
# Run Stata code first

# load data
    require(foreign)
    # survey data
    data <- read.dta("Female18-30.temp")
    attach(data)
    # data on UNMIL uniformed personnel
    unmil <- read.dta("UNMIL.dta")
    attach(unmil)

# Figure 1
    # compute proportions receiving different amounts from UN and non-UN transactional sex partners
    get.share <- function(val) {
        val <- table(val)
        lab <- names(val)
        names(val)[lab=="More than 200 dollars"] <- "$200+"
        names(val)[lab=="200 US dollars"] <- "$200"
        names(val)[lab=="100 US dollars"] <- "$100"
        names(val)[lab=="Fifty US dollars"] <- "$50"
        names(val)[lab=="Twenty US dollars"] <- "$20"
        names(val)[lab=="Ten US dollars"] <- "$10"
        names(val)[lab=="Five US dollars"] <- "$5"
        names(val)[lab=="One US dollar"] <- "$1"
        val <- val[lab!="Not applicable" & lab!="Don't know" & lab!="Refuse"]
        lab <- names(val)
        val <- c(val[lab=="$1"], val[lab=="$5"], val[lab=="$10"], val[lab=="$20"], val[lab=="$50"], val[lab=="$100"], val[lab=="$200"], val[lab=="$200+"])
        val <- val/sum(val)
    }
    un.plot <- get.share(tsvaluelast[tsUNlast=="Yes"])
    no.plot <- get.share(tsvaluelast[tsUNlast=="No"])
    # expand vectors for plotting
    no.plot <- c(no.plot[1], rep(0, 10), no.plot[2], rep(0, 10), no.plot[3], rep(0, 10), no.plot[4], rep(0, 10), no.plot[5], rep(0, 10), no.plot[6], rep(0, 10), no.plot[7], rep(0, 10), no.plot[8], rep(0, 3))
    un.plot <- c(0, 0, un.plot)
    un.plot <- c(rep(0, 2), un.plot[1], rep(0, 10), un.plot[2], rep(0, 10), un.plot[3], rep(0, 10), un.plot[4], rep(0, 10), un.plot[5], rep(0, 10), un.plot[6], rep(0, 10), un.plot[7], rep(0, 10), un.plot[8])
    label <- c(rep("", 2), "$1", rep("", 10), "$5", rep("", 10), "$10", rep("", 10), "$20", rep("", 10), "$50", rep("", 10), "$100", rep("", 10), "$200", rep("", 10), "$200+", rep("", 1))
    # draw bar plots
    x <- barplot(no.plot, ylim=c(0, .35), ylab="Share of most recent transactions", xlab="Value of anything received in most recent transaction (in USD)", space=-.75, col="grey40", names.arg="")
    barplot(c(0, 0, un.plot), add=T, space=-.75, col="grey70")
    axis(1, x[label!=""], label[label!=""])
    legend(14.5, .31, c("Man worked for UN", "Man did not work for UN"), fill=c("grey70", "grey40"))

# Figure 3
# bar plot of yearly share of women at risk who engage in transactional sex for the first time in that year

    # note that since women interviewed in 2012 were ages 18-30, age distribution of risk set varies substantially over time
    # need to make sure that women included in risk set for any given year represent the same age range
    # this age range for which we have data available depends on years we want to assess
    # for example, for 2002 to 2008, we have data for women aged 14-20 in a given year, because:
    # - women aged 20 in 2002 are 30 at the time of the survey in 2012, i.e. at the upper age limit for sensitive interviews
    # - women aged 14 in 2008 are 18 at the time of the survey in 2012, i.e. at the lower age limit for sensitive interviews

    # for bar plot, use the range of years that is optimal in terms of covering the largest total annual cases of women at risk, given that it:
    # - includes the year of UN deployment (2003) and the year right before and right after deployment,
    # - includes women aged 16-18 in the risk set, because they are particularly likely to engage in transactional sex for the first time
    # - does not include girls aged 13 or younger in the risk set, because they are relatively unlikely to engage in transactional sex

    # range of possible pre-deployment years to include
    pre <- 1996:2002
    # range of possible post-deployment years to include
    post <- 2004:2011
    # initialize grid computation
    totalcases <- c()
    out <- vector("list", length(post) * length(pre))
    count <- 1
    # for all possible ranges of years, compute age range for which we have data available, yearly size of risk set, and yearly share engaging in first transactional sex
    for(i in 1:length(pre)) {
        for(j in 1:length(post)) {
            .sampleatrisk <- .adjyearlyshare <- c()
            .fail <- F
            .years <- c(pre[i:length(pre)], 2003, post[1:j])
            # upper bound of age range for which data available
            .upper <- 30 - (2012 - min(.years))
            # lower bound of age range for which data available
            .lower <- 18 - (2012 - max(.years))
            # make sure age range meets requirements
            if(.lower<=13 | .lower>16 | .upper<18) .fail <- T
            for(k in .years) {
                .age <- age - (2012 - k)
                # make sure we did in fact observe data across age range
                if(all(.age!=.upper) | all(.age!=.lower)) .fail <- T
                .mark <- .age >= .lower & .age <= .upper
                # size of risk set
                .atrisk <- sum(.mark & (.age <= tsfirstnum | tsany==0), na.rm=T)
                # number that engaged in first transactional sex
                .new <- sum(.age==tsfirstnum, na.rm=T)
                .adjyearlyshare <- c(.adjyearlyshare, .new/.atrisk)
                .sampleatrisk <- c(.sampleatrisk, .atrisk)
            }
            if(.fail) {
                out[[count]] <- NA
                totalcases <- c(totalcases, NA)
                count <- count + 1
                next
            }
            .adjyearlyshare <- cbind(years=.years, share=.adjyearlyshare, sample=.sampleatrisk, lower=.lower, upper=.upper)
            out[[count]] <- .adjyearlyshare
            totalcases <- c(totalcases, sum(.sampleatrisk))
            count <- count + 1
        }
    }
    # range of years that includes largest total number of women at risk
    to.plot <- out[[which.max(totalcases)]]

    # draw bar plot
    par(mar=c(3, 4.1, 3, 4.1))
    yincrements <- .05
    ymax <- .35
    # to determine vertical axis automatically:
    # ymax <- max(yincrements * (1 + (to.plot[, "share"] %/% yincrements)))
    xvals <- barplot(to.plot[, "share"], names.arg=to.plot[, "years"], ylim=c(0, ymax), ylab=paste("Bars: Share of women aged", to.plot[1, "lower"], "to", to.plot[1, "upper"], "that engage in first transactional sex"))
    # add troop levels
    lines(xvals[to.plot[, "years"] %in% year], (ymax * total/max(total))[year %in% to.plot[, "years"]])
    lines(rep(xvals[4], 2), c(0, (ymax * total/max(total))[1]), lty=2)
    yunmil <- seq(0, 16000, 3200)
    axis(4, ymax * yunmil/max(total), yunmil)
    mtext("Line: UNMIL uniformed personnel", 4, 3)

